//-------------------------------------------------------------------------//
//                                                                         //
//  This benchmark is an OpenMP C version of the NPB SP code. This OpenMP  //
//  C version is developed by the Center for Manycore Programming at Seoul //
//  National University and derived from the OpenMP Fortran versions in    //
//  "NPB3.3-OMP" developed by NAS.                                         //
//                                                                         //
//  Permission to use, copy, distribute and modify this software for any   //
//  purpose with or without fee is hereby granted. This software is        //
//  provided "as is" without express or implied warranty.                  //
//                                                                         //
//  Information on NPB 3.3, including the technical report, the original   //
//  specifications, source code, results and information on how to submit  //
//  new results, is available at:                                          //
//                                                                         //
//           http://www.nas.nasa.gov/Software/NPB/                         //
//                                                                         //
//  Send comments or suggestions for this OpenMP C version to              //
//  cmp@aces.snu.ac.kr                                                     //
//                                                                         //
//          Center for Manycore Programming                                //
//          School of Computer Science and Engineering                     //
//          Seoul National University                                      //
//          Seoul 151-744, Korea                                           //
//                                                                         //
//          E-mail:  cmp@aces.snu.ac.kr                                    //
//                                                                         //
//-------------------------------------------------------------------------//

//-------------------------------------------------------------------------//
// Authors: Sangmin Seo, Jungwon Kim, Jun Lee, Jeongho Nah, Gangwon Jo,    //
//          and Jaejin Lee                                                 //
//-------------------------------------------------------------------------//

#include <math.h>
#include "header.h"
#include "mpinpb.h"


void compute_rhs()
{
  //printf("hello in rhs.c\n");
  int i, j, k, m, kstartbdy, kendbdy;
  double aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1, wijk, wp1, wm1;
  int scount, rcount,kstart, kend,kgap,jgap,jstart,jend, krem;
  int scount2[total_nodes], rcount2[total_nodes], displ2[total_nodes];
  int t1,t2,t3;
  MPI_Status status;
  MPI_Request sendrequest, recvrequest;
  if (timeron) timer_start(t_rhs);
  //#pragma omp parallel default(shared) private(i,j,k,m,rho_inv,aux,uijk, \
                                       up1,um1,vijk,vp1,vm1,wijk,wp1,wm1)
  {
  //---------------------------------------------------------------------
  // compute the reciprocal of density, and the kinetic energy, 
  // and the speed of sound. 
  //---------------------------------------------------------------------
  //#pragma omp for schedule(static) nowait
  kgap= grid_points[2]/ total_nodes;
  krem = grid_points[2]% total_nodes;
  if(krem == 0){
	kstart=node*kgap;
	kend= kstart+kgap;}
  else {
	  if(node+1 <= krem) 
	  {
		kstart = node*(kgap+1); //1 , 4
		kend= kstart+kgap+1; // 4, 7,
	  }
	  else
	  {
		kstart= krem*(kgap+1)+(node - krem)*kgap; //	7 , 9
		kend = kstart +kgap; // 9, 11
	  }
	}

  for(i=0; i< total_nodes; i++)
  {
    if(i+1 <= krem){
	rcount2[i] = (kgap+1)*(JMAXP+1)*(IMAXP+1);
	scount2[i]= rcount2[i];}	
    else{
	rcount2[i]= kgap*(JMAXP+1)*(IMAXP+1);
	scount2[i]= rcount2[i];}
   if(i == 0)
	displ2[i] = 0; 
   else 
	displ2[i] = displ2[i-1]+rcount2[i-1];
   	
  }
  
  for (k = kstart; k < kend ; k++) {
    for (j = 0; j <= grid_points[1]-1; j++) {
      for (i = 0; i <= grid_points[0]-1; i++) {
        rho_inv = 1.0/u[k][j][i][0];
        rho_i[k][j][i] = rho_inv;
        us[k][j][i] = u[k][j][i][1] * rho_inv;
        vs[k][j][i] = u[k][j][i][2] * rho_inv;
        ws[k][j][i] = u[k][j][i][3] * rho_inv;
        square[k][j][i] = 0.5* (
            u[k][j][i][1]*u[k][j][i][1] + 
            u[k][j][i][2]*u[k][j][i][2] +
            u[k][j][i][3]*u[k][j][i][3] ) * rho_inv;
        qs[k][j][i] = square[k][j][i] * rho_inv;
        //-------------------------------------------------------------------
        // (don't need speed and ainx until the lhs computation)
        //-------------------------------------------------------------------
        aux = c1c2*rho_inv* (u[k][j][i][4] - square[k][j][i]);
        speed[k][j][i] = sqrt(aux);
      }
    }
  }
  MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&rho_i[0][0][0],rcount2,displ2, MPI_DOUBLE,MPI_COMM_WORLD); 
  MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&ws[0][0][0],rcount2,displ2, MPI_DOUBLE,MPI_COMM_WORLD); 
 	
   //---------------------------------------------------------------------
  // copy the exact forcing term to the right hand side;  because 
  // this forcing term is known, we can store it on the whole grid
  // including the boundary                   
  //---------------------------------------------------------------------
  //#pragma omp for schedule(static)
  for (k = kstart; k < kend; k++) {
    for (j = 0; j <= ny2+1; j++) {
      for (i = 0; i <= nx2+1; i++) {
        for (m = 0; m < 5; m++) {
          rhs[k][j][i][m] = forcing[k][j][i][m];
        }
      }
    }
  }
 
//us array k+1
  if(node > 0) MPI_Isend(&us[kstart][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node-1, 1, MPI_COMM_WORLD, &sendrequest);

//us array k-1 
 if (node > 0)  MPI_Irecv(&us[kstart-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE,node-1,2, MPI_COMM_WORLD,&recvrequest);


  //---------------------------------------------------------------------
  // compute xi-direction fluxes 
  //---------------------------------------------------------------------
 // #pragma omp master
  if(node == root)
 	 if (timeron) timer_start(t_rhsx);
 if(node == root)
	kstart = kstart +1;
 if(node == total_nodes-1)
	kend = kend-1;
 //#pragma omp for schedule(static) nowait
 //t1= MPI_Wtime();
  for (k = kstart; k < kend; k++) {
    for (j = 1; j <= ny2; j++) {
      for (i = 1; i <= nx2; i++) {
        uijk = us[k][j][i];
        up1  = us[k][j][i+1];
        um1  = us[k][j][i-1];

        rhs[k][j][i][0] = rhs[k][j][i][0] + dx1tx1 * 
          (u[k][j][i+1][0] - 2.0*u[k][j][i][0] + u[k][j][i-1][0]) -
          tx2 * (u[k][j][i+1][1] - u[k][j][i-1][1]);

        rhs[k][j][i][1] = rhs[k][j][i][1] + dx2tx1 * 
          (u[k][j][i+1][1] - 2.0*u[k][j][i][1] + u[k][j][i-1][1]) +
          xxcon2*con43 * (up1 - 2.0*uijk + um1) -
          tx2 * (u[k][j][i+1][1]*up1 - u[k][j][i-1][1]*um1 +
                (u[k][j][i+1][4] - square[k][j][i+1] -
                 u[k][j][i-1][4] + square[k][j][i-1]) * c2);

        rhs[k][j][i][2] = rhs[k][j][i][2] + dx3tx1 * 
          (u[k][j][i+1][2] - 2.0*u[k][j][i][2] + u[k][j][i-1][2]) +
          xxcon2 * (vs[k][j][i+1] - 2.0*vs[k][j][i] + vs[k][j][i-1]) -
          tx2 * (u[k][j][i+1][2]*up1 - u[k][j][i-1][2]*um1);

        rhs[k][j][i][3] = rhs[k][j][i][3] + dx4tx1 * 
          (u[k][j][i+1][3] - 2.0*u[k][j][i][3] + u[k][j][i-1][3]) +
          xxcon2 * (ws[k][j][i+1] - 2.0*ws[k][j][i] + ws[k][j][i-1]) -
          tx2 * (u[k][j][i+1][3]*up1 - u[k][j][i-1][3]*um1);

        rhs[k][j][i][4] = rhs[k][j][i][4] + dx5tx1 * 
          (u[k][j][i+1][4] - 2.0*u[k][j][i][4] + u[k][j][i-1][4]) +
          xxcon3 * (qs[k][j][i+1] - 2.0*qs[k][j][i] + qs[k][j][i-1]) +
          xxcon4 * (up1*up1 -       2.0*uijk*uijk + um1*um1) +
          xxcon5 * (u[k][j][i+1][4]*rho_i[k][j][i+1] - 
                2.0*u[k][j][i][4]*rho_i[k][j][i] +
                    u[k][j][i-1][4]*rho_i[k][j][i-1]) -
          tx2 * ( (c1*u[k][j][i+1][4] - c2*square[k][j][i+1])*up1 -
                  (c1*u[k][j][i-1][4] - c2*square[k][j][i-1])*um1 );
      }
    }

    //---------------------------------------------------------------------
    // add fourth order xi-direction dissipation               
    //---------------------------------------------------------------------
    for (j = 1; j <= ny2; j++) {
      i = 1;
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m]- dssp * 
          (5.0*u[k][j][i][m] - 4.0*u[k][j][i+1][m] + u[k][j][i+2][m]);
      }

      i = 2;
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp * 
          (-4.0*u[k][j][i-1][m] + 6.0*u[k][j][i][m] -
            4.0*u[k][j][i+1][m] + u[k][j][i+2][m]);
      }
    }

    for (j = 1; j <= ny2; j++) {
      for (i = 3; i <= nx2-2; i++) {
        for (m = 0; m < 5; m++) {
          rhs[k][j][i][m] = rhs[k][j][i][m] - dssp * 
            ( u[k][j][i-2][m] - 4.0*u[k][j][i-1][m] + 
            6.0*u[k][j][i][m] - 4.0*u[k][j][i+1][m] + 
              u[k][j][i+2][m] );
        }
      }
    }

    for (j = 1; j <= ny2; j++) {
      i = nx2-1;
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp *
          ( u[k][j][i-2][m] - 4.0*u[k][j][i-1][m] + 
          6.0*u[k][j][i][m] - 4.0*u[k][j][i+1][m] );
      }

      i = nx2;
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp *
          ( u[k][j][i-2][m] - 4.0*u[k][j][i-1][m] + 5.0*u[k][j][i][m] );
      }
    }
  }
  //rloop1 =rloop1 + MPI_Wtime() - t1;
  if(node == root)
  {
  if (timeron) timer_stop(t_rhsx);

  //---------------------------------------------------------------------
  // compute eta-direction fluxes 
  //---------------------------------------------------------------------
  if (timeron) timer_start(t_rhsy);
  }

  //us array k+1 recv
  if(node < total_nodes -1) MPI_Irecv(&us[kend][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,1, MPI_COMM_WORLD,&recvrequest);
  //us array k-1 send
  if(node < total_nodes -1)  MPI_Isend(&us[kend-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,2, MPI_COMM_WORLD, &sendrequest);

//vs array k+1
  if(node > 0) MPI_Isend(&vs[kstart][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node-1, 3, MPI_COMM_WORLD, &sendrequest);
  if(node < total_nodes -1) MPI_Irecv(&vs[kend][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,3, MPI_COMM_WORLD,&recvrequest);

//vs array k-1
 if(node < total_nodes -1)  MPI_Isend(&vs[kend-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,4, MPI_COMM_WORLD,&sendrequest);
 if (node > 0)  MPI_Irecv(&vs[kstart-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE,node-1,4, MPI_COMM_WORLD,&recvrequest);

 //t2= MPI_Wtime();
  //#pragma omp for schedule(static)
  for (k = kstart; k < kend; k++) {
    for (j = 1; j <= ny2; j++) {
      for (i = 1; i <= nx2; i++) {
        vijk = vs[k][j][i];
        vp1  = vs[k][j+1][i];
        vm1  = vs[k][j-1][i];

        rhs[k][j][i][0] = rhs[k][j][i][0] + dy1ty1 * 
          (u[k][j+1][i][0] - 2.0*u[k][j][i][0] + u[k][j-1][i][0]) -
          ty2 * (u[k][j+1][i][2] - u[k][j-1][i][2]);

        rhs[k][j][i][1] = rhs[k][j][i][1] + dy2ty1 * 
          (u[k][j+1][i][1] - 2.0*u[k][j][i][1] + u[k][j-1][i][1]) +
          yycon2 * (us[k][j+1][i] - 2.0*us[k][j][i] + us[k][j-1][i]) -
          ty2 * (u[k][j+1][i][1]*vp1 - u[k][j-1][i][1]*vm1);

        rhs[k][j][i][2] = rhs[k][j][i][2] + dy3ty1 * 
          (u[k][j+1][i][2] - 2.0*u[k][j][i][2] + u[k][j-1][i][2]) +
          yycon2*con43 * (vp1 - 2.0*vijk + vm1) -
          ty2 * (u[k][j+1][i][2]*vp1 - u[k][j-1][i][2]*vm1 +
                (u[k][j+1][i][4] - square[k][j+1][i] - 
                 u[k][j-1][i][4] + square[k][j-1][i]) * c2);

        rhs[k][j][i][3] = rhs[k][j][i][3] + dy4ty1 * 
          (u[k][j+1][i][3] - 2.0*u[k][j][i][3] + u[k][j-1][i][3]) +
          yycon2 * (ws[k][j+1][i] - 2.0*ws[k][j][i] + ws[k][j-1][i]) -
          ty2 * (u[k][j+1][i][3]*vp1 - u[k][j-1][i][3]*vm1);

        rhs[k][j][i][4] = rhs[k][j][i][4] + dy5ty1 * 
          (u[k][j+1][i][4] - 2.0*u[k][j][i][4] + u[k][j-1][i][4]) +
          yycon3 * (qs[k][j+1][i] - 2.0*qs[k][j][i] + qs[k][j-1][i]) +
          yycon4 * (vp1*vp1       - 2.0*vijk*vijk + vm1*vm1) +
          yycon5 * (u[k][j+1][i][4]*rho_i[k][j+1][i] - 
                  2.0*u[k][j][i][4]*rho_i[k][j][i] +
                    u[k][j-1][i][4]*rho_i[k][j-1][i]) -
          ty2 * ((c1*u[k][j+1][i][4] - c2*square[k][j+1][i]) * vp1 -
                 (c1*u[k][j-1][i][4] - c2*square[k][j-1][i]) * vm1);
      }
    }

    //---------------------------------------------------------------------
    // add fourth order eta-direction dissipation         
    //---------------------------------------------------------------------
    j = 1;
    for (i = 1; i <= nx2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m]- dssp * 
          ( 5.0*u[k][j][i][m] - 4.0*u[k][j+1][i][m] + u[k][j+2][i][m]);
      }
    }

    j = 2;
    for (i = 1; i <= nx2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp * 
          (-4.0*u[k][j-1][i][m] + 6.0*u[k][j][i][m] -
            4.0*u[k][j+1][i][m] + u[k][j+2][i][m]);
      }
    }

    for (j = 3; j <= ny2-2; j++) {
      for (i = 1; i <= nx2; i++) {
        for (m = 0; m < 5; m++) {
          rhs[k][j][i][m] = rhs[k][j][i][m] - dssp * 
            ( u[k][j-2][i][m] - 4.0*u[k][j-1][i][m] + 
            6.0*u[k][j][i][m] - 4.0*u[k][j+1][i][m] + 
              u[k][j+2][i][m] );
        }
      }
    }

    j = ny2-1;
    for (i = 1; i <= nx2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp *
          ( u[k][j-2][i][m] - 4.0*u[k][j-1][i][m] + 
          6.0*u[k][j][i][m] - 4.0*u[k][j+1][i][m] );
      }
    }

    j = ny2;
    for (i = 1; i <= nx2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp *
          ( u[k][j-2][i][m] - 4.0*u[k][j-1][i][m] + 5.0*u[k][j][i][m] );
      }
    }
  }
  //rloop2 =rloop2 + MPI_Wtime() - t2;

  //#pragma omp master
  if(node == root)
  {
  if (timeron) timer_stop(t_rhsy);

  //---------------------------------------------------------------------
  // compute zeta-direction fluxes 
  //---------------------------------------------------------------------
  if (timeron) timer_start(t_rhsz);
  }



//qs array k+1
  if(node > 0) MPI_Isend(&qs[kstart][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node-1, 5, MPI_COMM_WORLD,&sendrequest);
  if(node < total_nodes -1) MPI_Irecv(&qs[kend][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,5, MPI_COMM_WORLD,&recvrequest);

//qs array k-1
 if(node < total_nodes -1)  MPI_Isend(&qs[kend-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,6, MPI_COMM_WORLD,&sendrequest);
 if (node > 0)  MPI_Irecv(&qs[kstart-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE,node-1,6, MPI_COMM_WORLD,&recvrequest);

//square array k+1
  if(node > 0) MPI_Send(&square[kstart][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node-1, 7, MPI_COMM_WORLD);
  if(node < total_nodes -1) MPI_Recv(&square[kend][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,7, MPI_COMM_WORLD,&status);

//square array k-1
 if(node < total_nodes -1)  MPI_Send(&square[kend-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE, node+1,8, MPI_COMM_WORLD);
 if (node > 0)  MPI_Recv(&square[kstart-1][0][0], 1*(JMAXP+1)*(IMAXP+1),MPI_DOUBLE,node-1,8, MPI_COMM_WORLD,&status);
 
 /*if(total_nodes > 1){
 MPI_Wait(&sendrequest,&status);
 MPI_Wait(&recvrequest,&status);}*/

  //#pragma omp for schedule(static)
  t3= MPI_Wtime();
  for (k = kstart; k < kend; k++) {
    for (j = 1; j <= grid_points[1]-2; j++) {
      for (i = 1; i <= grid_points[0]-2; i++) {
        wijk = ws[k][j][i];
        wp1  = ws[k+1][j][i];
        wm1  = ws[k-1][j][i];

        rhs[k][j][i][0] = rhs[k][j][i][0] + dz1tz1 * 
          (u[k+1][j][i][0] - 2.0*u[k][j][i][0] + u[k-1][j][i][0]) -
          tz2 * (u[k+1][j][i][3] - u[k-1][j][i][3]);

        rhs[k][j][i][1] = rhs[k][j][i][1] + dz2tz1 * 
          (u[k+1][j][i][1] - 2.0*u[k][j][i][1] + u[k-1][j][i][1]) +
          zzcon2 * (us[k+1][j][i] - 2.0*us[k][j][i] + us[k-1][j][i]) -
          tz2 * (u[k+1][j][i][1]*wp1 - u[k-1][j][i][1]*wm1);

        rhs[k][j][i][2] = rhs[k][j][i][2] + dz3tz1 * 
          (u[k+1][j][i][2] - 2.0*u[k][j][i][2] + u[k-1][j][i][2]) +
          zzcon2 * (vs[k+1][j][i] - 2.0*vs[k][j][i] + vs[k-1][j][i]) -
          tz2 * (u[k+1][j][i][2]*wp1 - u[k-1][j][i][2]*wm1);

        rhs[k][j][i][3] = rhs[k][j][i][3] + dz4tz1 * 
          (u[k+1][j][i][3] - 2.0*u[k][j][i][3] + u[k-1][j][i][3]) +
          zzcon2*con43 * (wp1 - 2.0*wijk + wm1) -
          tz2 * (u[k+1][j][i][3]*wp1 - u[k-1][j][i][3]*wm1 +
                (u[k+1][j][i][4] - square[k+1][j][i] - 
                 u[k-1][j][i][4] + square[k-1][j][i]) * c2);

        rhs[k][j][i][4] = rhs[k][j][i][4] + dz5tz1 * 
          (u[k+1][j][i][4] - 2.0*u[k][j][i][4] + u[k-1][j][i][4]) +
          zzcon3 * (qs[k+1][j][i] - 2.0*qs[k][j][i] + qs[k-1][j][i]) +
          zzcon4 * (wp1*wp1 - 2.0*wijk*wijk + wm1*wm1) +
          zzcon5 * (u[k+1][j][i][4]*rho_i[k+1][j][i] - 
                  2.0*u[k][j][i][4]*rho_i[k][j][i] +
                    u[k-1][j][i][4]*rho_i[k-1][j][i]) -
          tz2 * ((c1*u[k+1][j][i][4] - c2*square[k+1][j][i])*wp1 -
                 (c1*u[k-1][j][i][4] - c2*square[k-1][j][i])*wm1);
      }
    }
  }

  k = 1;
  if(node == (int) k/kgap){
  //#pragma omp for schedule(static) nowait
  for (j = 1; j <= grid_points[1]-2; j++) {
    for (i = 1; i <= grid_points[0]-2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m]- dssp * 
          (5.0*u[k][j][i][m] - 4.0*u[k+1][j][i][m] + u[k+2][j][i][m]);
      }
    }
  }
 }
  k = 2;
  if(node == (int) k/ kgap){
  //#pragma omp for schedule(static) nowait
  for (j = 1; j <= grid_points[1]-2; j++) {
    for (i = 1; i <= grid_points[0]-2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp * 
          (-4.0*u[k-1][j][i][m] + 6.0*u[k][j][i][m] -
            4.0*u[k+1][j][i][m] + u[k+2][j][i][m]);
      }
    }
  }
}
  
 kstartbdy = kstart;
 kendbdy= kend;
 if(kstartbdy < 3) kstartbdy =3;
 if(kendbdy > grid_points[2]-3) kendbdy = grid_points[2] -3;
 
  //#pragma omp for schedule(static) nowait
  //for (k = 3; k <= grid_points[2]-4; k++) {
   for (k = kstartbdy; k < kendbdy; k++) { 
    for (j = 1; j <= grid_points[1]-2; j++) {
      for (i = 1; i <= grid_points[0]-2; i++) {
        for (m = 0; m < 5; m++) {
          rhs[k][j][i][m] = rhs[k][j][i][m] - dssp * 
            ( u[k-2][j][i][m] - 4.0*u[k-1][j][i][m] + 
            6.0*u[k][j][i][m] - 4.0*u[k+1][j][i][m] + 
              u[k+2][j][i][m] );
        }
      }
    }
  }

  k = grid_points[2]-3;
 // #pragma omp for schedule(static) nowait
 if(node == (int) k / kgap){
  for (j = 1; j <= grid_points[1]-2; j++) {
    for (i = 1; i <= grid_points[0]-2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp *
          ( u[k-2][j][i][m] - 4.0*u[k-1][j][i][m] + 
          6.0*u[k][j][i][m] - 4.0*u[k+1][j][i][m] );
      }
    }
  }
}
 
 k = grid_points[2]-2;
  //#pragma omp for schedule(static)
 if(node == (int) k/ kgap){
  for (j = 1; j <= grid_points[1]-2; j++) {
    for (i = 1; i <= grid_points[0]-2; i++) {
      for (m = 0; m < 5; m++) {
        rhs[k][j][i][m] = rhs[k][j][i][m] - dssp *
          ( u[k-2][j][i][m] - 4.0*u[k-1][j][i][m] + 5.0*u[k][j][i][m] );
      }
    }
  }
}
  //#pragma omp master
  if(node ==root) 
	  if (timeron) timer_stop(t_rhsz);
 //rloop3 =rloop3 + MPI_Wtime() - t3;

  //#pragma omp for schedule(static) nowait
   for (k = kstart; k < kend; k++) {
    for (j = 1; j <= ny2; j++) {
      for (i = 1; i <= nx2; i++) {
        for (m = 0; m < 5; m++) {
          rhs[k][j][i][m] = rhs[k][j][i][m] * dt;
        }
      }
    }
  }
  
   MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&speed[0][0][0],rcount2,displ2, MPI_DOUBLE,MPI_COMM_WORLD); 

  } //end parallel
  if (timeron) timer_stop(t_rhs);
}
